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We present the first simulations of non-headon (grazing) 
collisions of binary black holes in which the black hole sin- 
gularities have been excised from the computational domain. 
Initially two equal mass black holes m are separated a dis- 
tance « 10m and with impact parameter « 2m. Initial data 
are based on superposed, boosted (velocity ~ 0.5c) solutions 
of single black holes in Kerr-Schild coordinates. Both rotating 
and non-rotating black holes are considered. The excised re- 
gions containing the singularities are specified by following the 
dynamics of apparent horizons. Evolutions of up to t ~ 35m 
are obtained in which two initially separate apparent horizons 
are present for t ~ 3.8m. At that time a single enveloping ap- 
parent horizon forms, indicating that the holes have merged. 
Apparent horizon area estimates suggest gravitational radia- 
tion of about 2.6% of the total mass. The evolutions end after 
a moderate amount of time because of instabilities. 

Introduction: Gravitational wave detectors |1| will 
soon begin searching for gravitational radiation from as- 
trophysical binary compact objects. To understand these 
observations, and to predict parameter regimes in which 
to search for their radiation, efforts are underway to 
model the interaction of compact sources. We report 
here a direct numerical simulation of interacting spin- 
ning black hole binaries, in genuinely hyperbolic (non- 
headon) trajectories. The initial spin angular momenta 
evolved here are either zero, or parallel to each other and 
perpendicular to the orbital plane. The interior of the 
equal mass holes and their interior singularities are ex- 
cised from the computation. (Our method is neither re- 
stricted to equal masses nor to parallel spins) . Evolution 
is carried out in a Cauchy scheme, in which the state of 
the gravitational system (the 3-spatial metric gab) and its 
rate of change (the 3-spatial extrinsic curvature Kab) are 
specified at one instant (i.e. on a 3-dimensional space- 
like hypersurface) and are then stepped to the next in- 
stant using an "ADM" |2] form of the Einstein evolution 
equations H. The evolution is unconstrained, and main- 
tenance of the constraint functions with small error is 
verified throughout the run. 

This work extends previous work on headon encoun- 
ters 0-0] • It is comparable to recent results of Briigmann 
pi: non-headon black hole evolution through to signifi- 
cant interaction and merger. But our approach has a 
novel feature: the singularity- excising character of the 
computation of generic encounters which allows "natu- 
ral" motion of the black holes through the computational 



grid. Singularity excision may be crucial to carrying out 
long term simulations predicting gravitational waveforms 
through several wave-cycles. 

Initial Data: We carry out three binary black hole simu- 
lations. Data is created with spinning holes, each of mass 
m, located at (±5m, ztm, 0), each with Kerr spin param- 
eter a. The holes are boosted in opposite x— directions 
with speed c/2, representing a grazing collision with im- 
pact parameter of 2m (and resulting total orbital angu- 
lar momentum in the z— direction). We distinguish three 
cases: case (I)- both holes have a = 0.5m opposite to the 
orbital angular momentum; case (II)- nonspinning holes 
a — 0; case (III)- both holes have a = 0.5m aligned with 
the total angular momentum. 

The total initial ADM mass of each simulation is 
2.31?7i, which agrees very well with the estimate given 
by the special relativistic limit ruADM = 27m, with 
7 = (1 - .52)-i/2 = 1.155, The total initial ADM an- 
gular momentum J = Jz is 0.0, l.llm?, and 2.34m^ for 
cases I, II, and III respectively (see |^]). 

The data setting technique is based on the boost- 
invariant Kerr-Schild [Q form of the Kerr black hole 
metric. Our Cauchy formulation requires first the solu- 
tion of the initial data problem. As outlined in [^l|-|3[ , 
superposed boosted Kerr-Schild data for two single holes 
produce a conformal background space; the physical data 
are solved via a York-conformal approach (solving four 
coupled elliptic equations) Q on this background. Note 
that even when an exact solution of the elliptic equa- 
tions is known, the error in the evolved solution will be 
determined by the inherent evolution-equation trunca- 
tion error. Therefore, the accuracy of elliptic solver em- 
ployed need just be consistent with this truncation error. 
For the discretization used here (Aa; = m/4) the trun- 
cation error is of order 5%. The quality of the data is 
validated by computing the constraints, normalized to 
a dimensionless quantity by the factor m^^ . Analyti- 
cally the constraints should be zero everywhere. In fact 
with the parameters of the problem, and with the cur- 
rent discretization and truncation error, the superposed 
background solution is acceptable with no further ellip- 
tic problem solution |1^ (i.e. the 0th order of the elliptic 
solver). However, as we progress to larger and better 
resolved evolutions, we will find it mandatory to cycle 
through the elliptic solve step ||l^ to obtain satisfactory 
solution of the constraints. Figure 1 presents the Hamil- 



tonian constraint for case III, evaluated at integration 
timestep t = Sm along the x— axis, together with a time 
history of the I2 norm (over volume outside the horizons, 
and excluding the outer boundary region) of the Hamil- 
tonian constraint and the similarly normed momentum 
constraint. The late time rise in the momentum con- 
straint in Figure 1 shows the beginning of the exponen- 
tial mode that appears at about t = 36m and ends the 
simulation. We have quite good constraint behavior, of 
order 0.4%, with peak errors in the Hamiltonian of order 
5% until that time. 

Evolution Methods: The time-evolutions presented 
here are done using AGAVE, a code that solves the Ein- 
stein equations in an ADM 3-1-1 form via finite difference 
techniques |^. A parallel implementation is obtained 
with the use of MPI jl^l, employing the Cactus com- 
putational toolkit |l7| solely to aid in this task. AGAVE 
is a major revision of the Binary Black Hole Grand Chal- 
lenge Alliance Cauchy code fljjl^ . The lapse function a 
and shift vector /?* express coordinate conditions which 
are chosen to allow the black holes to move freely. For 
our simulations, prior to the time that a single black 
hole surrounding the incoming pair is detected, we use 
a superposition of functions from boosted black holes: 
a = ai + a2 — 1 , /3' = /3| -I- /3| , where these functions are 
centered with the current location of the holes, and with 
the velocity initially obtained from Newtonian approxi- 
mation to the trajectories of the holes and subsequently 
inferred from the history of the locations of the appar- 
ent horizons(see below); after the detected merger, we 
use the lapse and shift of a single black hole with a mass 
which is the sum of the original bare masses, and angular 
momentum which is the (naive vectorial) sum of the spin 
and orbital angular momentum in the original system. 
(See Discussion, below.) 

The interior of the black holes is excised (Unruh, 
quoted in |Q). We use the apparent horizon surface, 
locatable at each time-slice, as a marker for the excision. 
We utilize a combination of two different finite difference 
methods to find the apparent horizon: a direct solver pi| , 
and a curvature flow method ||22| . Once the apparent 
horizon is located, we define a mask function that delin- 
eates the excluded region (interior to the holes) from the 
computation. The result is that we literally evolve two 
holes moving freely through the computational domain. 
That domain is a 161'^ lattice, corresponding at our reso- 
lution to a cube {AOm)^ (±20m in each direction from the 
centered origin). However, boundary conditions are set 
by providing Dirichlet boundary conditions for gat and 
blending | |23| , p^ outwards from a sphere of radius 19m the 
computational solution of Kab to an analytically given 
(time-dependent) solution for Kat at the outer boundary 
sphere. "Blending" means taking a linear combination 
of values from the computed and the analytically given 
solution, over a few (here, four) spatial zones, reducing 
gradients and second derivatives at the boundary. The 



analytic blending solution is created by superposition of 
boosted holes given by the initial data construction (with 
centers and velocities propagated according to the lapse 
and shift computation), or after the merger by the final 
estimated black hole with post merger lapse and shift. 

The discretization of the Einstein equations is consis- 
tent to second order accuracy. On the time scale where 
instabilities do not play a significant role, the convergence 
rate of this code is « 1.6, reduced from 2 apparently be- 
cause of extrapolation at the excision boundaries. 
Results: To the current accuracy of the code, cases I-HI 
behave similarly. The total proper area of the apparent 
horizon A for case (I) is shown in Figure 2. The value of 
A is particularly interesting since it provides a measure of 
the total mass contained in the apparent horizon. For a 
given black hole of mass m and spin parameter a its area 
is Abh ~ 47r(i?^-|-a^) (with i?+ — m+^/m? — a^). Since 
at early times there is no common apparent horizon the 
total area is approximately A « Abhi +Abh2 — "^Abhi, 
as the holes merge the total mass enclosed in the common 
horizon is (roughly) expected to double, and hence its 
area would be four times as bigger, ie. for a non-spining 
final black hole A w ATT{2{2m)f < AAbhi- Therefore, a 
plot of A vs. time (like the one in figure 2) shows a con- 
siderable 'jump' at the time the holes merge t « 3.8?n. 
Additionally, effects of the outer boundary can be clearly 
seen in figure 2. For a ±10 grid an abrupt 'kink' is seen 
aX t ~ lOm while in the ±20 grid the 'kink' appears at 
t « 20r7T,. At about t « 36r7i [t « 26m) apparent insta- 
bilities in the ±20m (±10m,) grid cause a rapid increase 
in the computed horizon size and eventually crash the 
run. Thus at t sa 35m the solution becomes untrustwor- 
thy. While the simulation is free of boundary effects the 
coincidence of the measured horizon area values supports 
confidence in the results. Figures 3A - 3F track the ap- 
parent horizons through the merger for case I. A single 
enveloping black hole appears at i « 3.8m,. The horizon 
oscillates and grows slightly. 

We have in place Cauchy-characteristic extraction, 
where the Cauchy solution sampled at some "large" ra- 
dius acts as data for a characteristic evolution to infinity 
p5| , p6| for waveform extraction. We also can compute 
the Newman Penrose tensor -04, which captures at null 
infinity the outgoing radiation. Additionally, we are de- 
veloping a perturbative radiation extraction module. We 
are preparing an article explaining how these tools are 
applied and illustrating the radiation patterns obtained 
from these simulations. 

Discussion and Future Directions: The simulations 
reported here are genuinely, but not excessively, hyper- 
bolic encounters. A Newtonian estimate gives a free fall 
velocity of 0.4c from infinity, as compared with the ve- 
locity 0.5c specified in our initial data. Future work will 
concentrate on generic hyperbolic and elliptic orbits. 

Ongoing research concerns the late-time stability of the 
black hole simulations. We have carried out a number of 



l-dimensional simulations, all of which have longer term 
stability than this 3-climensional simulation of merged 
holes. We are investigating the behavior of the differ- 
encing scheme at the inner boundary. (The one we use 
behaves well in the spherical case.) We are implementing 
a new outer boundary algorithm which has been shown to 
be robustly stable in a linearized version of the code p^ . 
We are developing more sophisticated gauges based on 
elliptic equations for the lapse and the shift. These in- 
clude the minimal distortion and minimal shear gauges 
p8[ , and other elliptic gauges ||ll|,^ . Stable evolution of 
single black holes is quite sensitive to gauge conditions, 
and we anticipate much useful science from future im- 
provement in the lifetime of our simulations of black hole 
mergers. 

Our gauge and boundary conditions for the final 
merged black hole naively assume that all the initial mass 
(i.e. M final — 2m) and angular momentum resides in the 
final hole: J final = a final x M final- For cases I, II, III 
our gauge takes afinai = (0,0.25,0.5) x M final- These 
estimates do not take into account the emission of en- 
ergy and angular momentum during the dynamics, nor 
the 7 factor in the initial mass and angular momentum. 
The actual post-collision mass and angular momentum of 
the residual hole will be evaluated to further improve the 
simulations; behavior of the code is robust under changes 
in the final assumed mass and spin. 

Of extreme interest is the size of the final apparent 
horizon. The total initial ADM mass leads to horizon 
area of 47r(2 x 2.31m)^ « 268m^. The post-merger nu- 
merically computed apparent horizon area (Figure 2) is 
about 255to^, 5% smaller than this estimate. This mea- 
sure would give a preliminary indication that total en- 
ergy radiated in this simulation is about 2.6%. However, 
we have yet to complete a 3-dimensional event horizon 
tracker, which will allow a correct comparison of the ini- 
tial and final event horizon area. 

The present work demonstrates the first simulation of 
binary black hole systems via the excision of singularities. 
The datasets evolved are not only useful for validation of 
the techniques employed here but as valid datasets in an 
astrophysical sense for the final "plunge" of the merger. 
In this work we: (a) demonstrate well behaved (conver- 
gent) descriptions of the black holes as they evolve; (b) 
show that apparent horizon tracking and black hole exci- 
sion can produce dynamical multi-black hole spacetimes, 
with reasonably well controlled errors for a considerable 
length of time (long enough for an accurate modeling 
of the merger phase); and (c) demonstrate that rela- 
tively unsophisticated gauge functions a and j3 can lead 
to physically interesting evolution lifetimes. 
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FIG. 1. For case I (and grid±20), the Hamiltonian and mo- 
mentum constraints, on the domain of outer communication 
(outside the apparent horizon(s) and inside the outer bound- 
ary blending zone). We give the time history of the I2 (rms) 
norm of the Hamiltonian (solid line) and the I2 norm over all 
three components of the momentum constraint (dotted line) . 
The momentum I2 is constructed only along coordinate lines 
(all that is available from this computation); the Hamiltonian 
I2 is computed from the whole volume. The sudden change in 
the errors at t fa 4m occurs when a single outer apparent hori- 
zon envelops the merging holes. Also, the drop at t « 20m 
is due to boundary effects. The inset shows the Hamiltonian 
constraint along the a;— axis at time t = 3m.. 
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FIG. 2. The area of the apparent horizon(s) (transition to 
a single horizon at t ~ 4m) for case I. For a smaller domain 
(zblOm, dashed line) the simulation runs to i ~ 26m and 
exhibits strong boundary effects at i ~ 10m. In the larger 
(±20m) domain (solid line) boundary effects show at later 
time, around 20m,. Instabilities cause the measured area to 
rise abruptly at f « 36771 and eventually stop the simulation. 




FIG. 3. For case I, time history of the horizons. The times 
corresponding to figures 3A-3F are t = 0, 2.6m, 5.1m,, 8.8m, 
13.8m, 18.8m. These are coordinate plots; the correspond- 
ing areas appear in Figure 2. After the merger the horizon 
oscillates through a fraction of a cycle. 



